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Abstract 

A new method of post-Newtonian approximation (PNA) for weak 
gravitational fields is presented together with its application to test an 
alternative, scalar theory of gravitation. The new method consists in 
defining a one-parameter family of systems, by applying a Newtonian 
similarity transformation to the initial data that defines the system 
of interest. This method is rigorous. Its difference with the standard 
PNA is emphasized. In particular, the new method predicts that the 
internal structure of the bodies does have an influence on the motion 
of the mass centers. The translational equations of motion obtained 
with this method in the scalar theory are adjusted in the solar system, 
and compared with an ephemeris based on the standard PNA of GR. 



1 This is a revised and expanded version of the text of a lecture given at the 8 th Conf. 
"Physical Interpretations of Relativity Theory" (London, September 2002), organized by 
the British Society for the Philosophy of Sciences. The initial version will appear in the 
Proceedings of that Conference (M. C. Duffy, ed.). 
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1 Introduction 



One first expects from a theory of gravity that it should provide an ac- 
curate celestial mechanics. In other words, the theory should tell us how 
massive celestial bodies precisely move with respect to each other under the 
effect of the gravitational field produced by them all. Thus, Einstein's gen- 
eral relativity (GR) won its first advantage over the older theory of Newton 
when it gave an explanation to Mercury's residual advance in perihelion. In 
1972, Weinberg stated about this explanation ([TJ, p. 198): "This is by far the 
most important experimental verification of general relativity." It is hence 
extremely important for a theory of gravitation, not only that it produces 
accurate ephemerides, but even more that one is sure that it really produces 
those ephemerides, i.e., that the solution of the approximate equations used 
in the computation does approach accurately enough the relevant solution 
of the exact equations. The works of Fock |2j and Chandrasekhar [3j aimed 
at answering the latter question for GR. The later work on celestial mechan- 
ics in GR relies on essentially the same approximation scheme as these two 
works, which are equivalent in this regard. Yet in 1966, Synge, who did 
know Fock's work (which is quoted in Ref. jl]) and most probably knew also 
that of Chandrasekhar, wrote 0:"! am still waiting for a rational treatment 
of the dynamics of the solar system according to Einstein's theory. In the 
very nature of the case, any argument must be of an approximate nature; an 
assessment of the error is a primary desideratum." Comparing his successive 
sentences, we may infer that Synge was not satisfied with the approximation 
method used in the works El nor with the one he himself proposed with 
coworkers [HI E] , and which was limited to stationary fields — this restriction 
is indeed inappropriate to describe the solar system in a realistic way. 

The aim of this paper is to summarize the principles, the development 
and the numerical implementation of a new approximation method for ce- 
lestial mechanics in relativistic theories of gravitation. This approximation 
method might have satisfied Synge, perhaps, because it is mathematically 
sound and general, and because it too predicts a salient result which he 
found in his work for GR |H1 E], namely the fact that, in such theories, the 
internal structure of a body does influence the gravitational field produced 
by it, hence also the motion of external bodies CHI- The new method 
consists basically in associating a one-parameter family (S>,) of gravitating 
systems with the physically given system S, by defining a family of initial 
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conditions. It was initiated by Futamase & Schutz [7j for GR, with fur- 
ther mathematical developments given by Rendall jH]- However, Futamase 
& Schutz |7| assumed a very restrictive initial condition for the spatial metric. 
As to Rendall jH], he considered an a priori given one-parameter family of 
solutions of the field equations, without investigating the definition of a such 
family from the given system S. Moreover, these two works were limited to 
the local equations and some of their mathematical properties. In particular, 
they did not provide equations of motion for the mass centers of a system 
of extended bodies, as one needs to compute an ephemeris. We came to the 
new method independently [HIE], to test an alternative theory of gravita- 
tion based on just a scalar field [El H3j, and we did obtain such equations of 
motion ^3 El CH|. That scalar theory gives the same predictions as GR for 
light rays [TI]. Therefore, it is worth testing this theory further. Moreover, 
since that theory is much simpler than GR, it is easier to implement the new 
method for that theory, as well as to discuss the difference between the new 
method and the standard PNA. 



2 General framework: the method of asymp- 
totic expansions 

As is well-known, an asymptotic expansion of a real function (p of the real 
variable A in the neighborhood of some value A is an expression 

(p{\) = a ^o(A) + ... + a„^(A) + R{\), (1) 

the known functions ifi®, ip n being positive and belonging to a definite com- 
parison set E, endowed with certain properties, and with ip ^> ... tp n R 
as A — > A ^UJ. In physics, the relevant value is usually A = 0, and one speaks 
of the "small parameter" A. In particular, if the behaviour as A — > is regu- 
lar enough, a Taylor expansion may apply, so that ^(A) = \ k {k = 0, ...,n). 
However, it may be that the Taylor expansion can be pushed only to some 
order n = n max , beyond which a more accurate expansion can be obtained 
only if one accepts to consider more general functions, e.g. ones involving 
a fractional exponent k. This remark is relevant to weak-field expansions in 
relativistic theories of gravitation. 

Now consider a boundary-value problem defined for a given system of 
partial differential equations (PDE's), and assume that a small parameter 
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A can be denned for this problem, which means in fact that a family (Pa) 
of problems can be defined. The method of asymptotic expansions for this 
problem consists in trying to write each scalar component of the solution 
of Pa, say tp^ (where i — 1, ...,m with m the number of scalar unknowns 
involved in the system of PDE's), as an asymptotic expansion in A. That 
expansion is thus assumed valid at each given point IgD, where D is the 
relevant domain for the independent variables (space and time, say), which 
are collectively denoted by X. Actually this domain itself may well depend on 
A, but, in order that one may write definite expansions (0) involving known 
functions ^fc(A), it is preferable to absorb this dependence in a redefinition of 
the independent variables such that the expansions indeed apply to any given 
point X in a domain independent of A. If we look for a Taylor expansion, we 
write thus: 

VX e D, ipf\x) = <Pm(X) + ... + y in {X)\ n + Ri(X, A), (2) 

with Ri(X, A) <C A™ as A — > 0. The problem which is really of physical 
interest, P, e.g. the initial-value problem for some (assumed) isolated self- 
gravitating system, is assumed to correspond to a given, small value Ao of the 
parameter. One of the difficulties of the method is definition of an adequate 
family (Pa), from the given problem P. Once this has been done and once 
expansions as A — > 0, and corresponding expanded equations, have been ob- 
tained for the family (Pa), they are then used for the finite value Ao- This 
means that the error involved is the value for Ao of the unknown remainder 
Ri(X,X). In the case of a Taylor expansion, however, the remainder will 
usually be 0(A n+1 ), and this uniformly with respect to X e D (the relevant 
domain D being often compact). Although this is only an asymptotic error 
estimate, thus not a numerical one, it can be said that, if A n+1 is negligible 
with respect to the experimental accuracy, then the n-th order expansion (J2J) 
is very likely to be enough accurate for a meaningful experimental test. This 
seems to be the best that one can hope in the current state of the theory of 
PDE's. 

Why is it useful to write expansions at all? Mainly because the resulting 
equations are very greatly simplified: as it turns out, all non-linearities are 
reported in the equations of the order zero, and those are often much simpler 
than the starting equations. We emphasize two points: 
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i) It is indeed necessary to define a family of boundary-value problems 
(instead of contenting oneself with just that problem which one is interested 
in). This is in order that it just make sense to try expansions in A: if we 
have defined a such family, we can then expand with respect to A each of 
the various coefficients that define the system of PDE's, and we also can 
expand the boundary values. Only in that case can we derive expansions of 
the solution fields, and expanded equations for them, and then solve these 
equations using the expanded boundary conditions. 

ii) One should indeed expand all independent unknown fields, (p\ for 
i = 1, ...,m (and not just those which one likes to expand). For this allows 
one to write each equation as a sum of terms, each of a definite order in 
A, this allowing in turn to separate the equations of the different orders: 
k — 0, ...,n — which is necessary, because an expansion like (j2J) means that a 
field ifl is, after expansion, split into the (n + 1) fields ip i0 , ...,ip in . Thus, 
if we write expansions for the m independent unknown fields, using expan- 
sions with (n + 1) terms (e.g. Taylor expansions of order n), then we get 
each of the m independent equations split into (n + 1) equations, so that 
we now have m(n + 1) equations for m(n + 1) unknowns. Whereas, if one 
expands only p among the independent fields (with p < m), using expan- 
sions with (n + 1) equations, then one has (of course!) no reason to split 
the equations. But if one nevertheless would do so, then one would have too 
much equations (if n > 1): m(n + l) equations for p(n + 1) + m — p unknowns. 

The above-described method is merely the general formulation of a nat- 
ural perturbation method for a system of PDE's, and it is certainly not new. 
Of course there is a vast literature on perturbation methods (see e.g. Refs. 
17-18 and references therein), and there are many common points between 
the different approaches. Yet we have not been able to find a description 
like the foregoing one, which applies very closely to what we actually did 
for the scalar theory of gravitation. Certainly also, the knowledge of PDE's 
and of rigorous perturbation methods has considerably improved since the 
fifties. This may explain why the method developed for weak fields in GR 
by Fock j2] and Chandrasekhar [3] does not fulfil the requirements i) and ii) 
above — in fact it is based on formally taking 1/c 2 as a small parameter (with 
c the velocity of light), and in expanding the gravitational field, but not the 
matter fields, in powers of 1/c 2 . 
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3 Asymptotic post-Newtonian approximation 
of the scalar theory 

It is the application of the foregoing method to the case of that "relativis- 
tic" theory of gravitation. 2 In such theories, the relevant boundary-value 
problem is the initial-value problem. This is due to the hyperbolic character 
of the gravitational equation, in other words it comes from the fact that, in 
such theories, gravitation propagates with a finite velocity, usually equal to 
the velocity of light c or close to c. This applies |T2| to the scalar theory 
investigated by the author. Note, however, that the complete system of local 
equations is not closed, hence in particular cannot be qualified "hyperbolic" , 
until one has postulated a definite behaviour for matter, by assuming a con- 
stitutive equation giving the material energy-momentum tensor T in terms of 
some matter fields. For the sake of simplicity, we assume a perfect barotropic 
fluid, for which one has only the pressure p and the (spatial) velocity u as 
independent matter fields. 3 

Thus we consider a given, isolated gravitating system S, made of N sepa- 
rated bodies. That system is defined by the barotropic state equations in the 
different bodies, and by the initial data for the independent matter fields p 
and u, as well as for the scalar gravitational field / and for its time derivative 
drf- (In GR, an initial data for the gravitational field is much more com- 
plicated, because the latter field is a tensor one, also because the Einstein 
equations are underdetermined, and above all because the initial conditions 
have to verify nonlinear constraint equations PT) I22|.) The first task is to de- 
fine a family (S\) of gravitating systems, by defining a family of initial data 
and state equations. We want to describe weakly gravitating systems, for 

2 The scalar theory investigated is indeed relativistic in the sense that it accounts for 
special relativity, and reduces to it if the gravitational constant G vanishes — but it is a 
preferred- frame theory. The summary of the scalar theory that is given in Ref. |12| . Sect. 
2, is sufficient (and not even necessary, we believe) for the present purpose. Being based 
on a scalar field, that theory is, of course, very different from the relativistic theory of 
gravitation proposed by Logunov et al. I2()j . However, both theories consider a flat 
"background" metric and a curved "effective" metric. 

3 The nature of the relevant boundary-value problem might a priori be expected to 
depend on the constitutive equation assumed. However, even if one assumed a very general 
matter behaviour, including anelasticity and dissipation effects, it seems that the initial- 
value problem would remain the "good problem" . The reason to believe this is that it is 
so for the heat equation, although the latter is parabolic. 
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which Newton's theory with its Euclidean space and absolute time must be 
an excellent approximation. For this to be true, the family (S\) must satisfy 
two conditions as A — > 0: i) the (physical, or "effective") space-time metric 
■y must tend towards a flat metric 7 , and ii) all fields must become closer 
and closer to "corresponding" Newtonian fields. For the scalar theory, which 
is bimetric, the flat metric 7 always coexists with the physical one 7, and 
their local difference is an increasing function of 1 — /, which is non-negative 
|12j . Condition i) is hence easy to be made precise for the scalar theory: 
it means simply that the maximum value of the field (1 — f^') must tend 
towards zero as A — > 0. The parameter is 

A = Su PxeM [l-/W(x,T = 0)]/2. (3) 

where M is the "space" manifold, i.e. the set of the positions x in the pre- 
ferred reference frame (PRF). 

To make condition ii) precise, we use a crucial result [12], namely 
the existence of an exact similarity transformation in Newton's theory for 
barotropic fluids. Suppose we have the following fields: pressure p^\ density 
= F^(p^), Newtonian potential uffl, and velocity u^, obeying the 
continuity equation, Poisson's equation, and Euler's equation. Then, for any 
A > 0, the fields pW(x,T) = AV x )(x, VX T), p (A) (x,T) = Ap (1) (x, \/\ T), 
f/^ A) (x,T) = Af4 1} (x, T), and U W(x,T) = uM(x, T), also obey 
these equations — provided the state equation for system Sa is F^ x \p^) = 
XF^ 1 ' (X~ 2 p^) . This similarity transformation defines the weak-field limit in 
Newton's theory itself: as A — > 0, the potential and the density in the bodies 
decrease like A (while the bodies keep the same size), the (orbital) velocities 
decrease like vA, and accordingly the time scale increases like 1/VA. 

This exact similarity transformation does not extend to a "relativistic" 
theory like the scalar theory, simply because the equations are different. In 
particular, the metric changes as the gravitational field becomes stronger, 
and this changes the motion in a complex manner. Recall, however, that 
we merely seek initial conditions for the fields in the scalar theory. It is 
then obvious that, in order to precise condition ii), we simply have to apply 
the Newtonian similarity transformation to the initial fields. For the matter 
fields, this application is indeed immediate, because for a perfect fluid these 
are common to Newton's theory and to a relativistic theory (via some slight 



7 



modifications): pressure p, proper rest-mass density p* (instead of the invari- 
ant density of Newton's theory), mass density of elastic energy II, coordinate 
velocity u = dx/dT (for the scalar theory, T is the preferred time coordinate 
or "absolute time", and x is the position in the preferred reference frame 
(PRF)). The mere difficulty is that, to apply the transformation, we must 
also associate a kind of "Newtonian potential" with the scalar gravitational 
field /. To do that, we assume in advance that a family (S\) of systems 
has been built, which is such that the orders in A of the matter fields are the 
same, as A — > 0, as in the Newtonian similarity transformation, and such that 
the field f( x ' admits an expansion of the form /( A ) = 1 + 0A fc + 0(\ k+1 ). We 
find then that one must have k = 1, i.e. the field 1 — is like A as A — > 0, 
thus like the potential in the Newtonian weak-field limit. 4 It is hence natural 
to define that field, more precisely V {x) = (c 2 /2)(l - /W), as the equiva- 
lent of the Newtonian potential uffl; the coefficient is needed to ensure that 
yW ~ as A — > 0. Moreover, the local difference between the flat metric 
7° and the curved physical metric 7 is an increasing function of 1 — /^(x, T). 

In that way, we get the initial conditions for system S\. Then we state 
powers expansions in A for the independent fields f,p, u, and we deduce 
expansions for the other fields. The set of the obtained expansions and ex- 
panded equations is internally consistent, moreover it is consistent with the 
initial aim, in that the equations of the order are indeed the "Euler-Newton" 
equations, i.e. the equations for a perfect fluid in Newton's theory. The first 
PN corrections correspond to the equations of the order 1. The details can 
be found in Ref. 12 (see Ref. 15 for a synopsis and a few complementary 
points). Let us mention here some important points: 

- We use a change of the mass and time units for system S\, adopting 
[M] A = A[M] and [T] A = [T]/\/A as the new units (where [M] and [T] are 
the starting units); then all fields are ord(A°), and the small parameter A is 
proportional to 1/c 2 (in fact A = (c /c) 2 , where c is the velocity of light in 
the starting units). That 1/c 2 , not 1/c, turns out to be the effective small 
parameter, is due to the fact that it is only 1/c 2 that enters the equations. 
Thus, the derivation of the expansions and expanded equations is straight- 
forward. In the standard PN scheme 1/c 2 is formally considered as a 

4 The small parameter considered in the present paper is A = e 2 , where e is that used 
in Ref. 12. 
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small parameter, and the matter fields are not expanded. 

- Recall, however, that such Taylor expansions in A (or 1/c 2 ) cannot be 
expected to apply to an arbitrary order n [see after Eq. (JTJ]. In GR, it has 
been shown by Rendall [H] that the attainable order n depends on the gauge 
condition: with the classical "harmonic gauge" , even the second PN approx- 
imation (n = 2) cannot be solved for a general spatially compact system 
of bodies — but a different gauge condition allows to reach the order n = 3. 
This kind of problem cannot occur in the scalar theory, because there is no 
gauge condition. However, only the first PN approximation (n = 1) has been 
studied in the scalar theory. If the matter sources are considered given (as 
was done in GR by Rendall |8j, then the PN gravitational field depends on 
integrals whose integrand is regular and vanishes outside the bodies [see Eqs. 
(f2"H) - (j2T)j) below], so the first PN approximation is unproblematic. 

- The equations of the order 1 are linear with respect to the fields of 
the order 1, so that the nonlinearities are confined to the zero-order (Euler- 
Newton) equations. The separation between the matter fields of the different 
orders, as well as the linearity in the PN fields, are characteristic of the 
asymptotic method, as contrasted with the standard PN scheme (21 EH- It 
seems worth to illustrate the difference between the two approximations by 
exhibiting some simple equations. 

4 Illustrating and commenting the difference 
between "standard" and "asymptotic" PNA's 

In the asymptotic PNA of the scalar theory, the rest-mass density in the 
PRF, pexact, is approximated as 5 

Pcxact = P(l)[l + 0(A 2 )], p(i) = p + Pl/c 2 , (4) 

and the like for the other fields, e.g. the velocity field u exa ct- (For the zero- 
order, Newtonian fields, we shall omit the index and shall thus keep the 
usual notation, while denoting henceforth the exact fields by the subscript 

5 The asymptotic expansion has the form QJi in the "varying units" defined above, 
but it has a slightly different form in invariable units |12| . 
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"exact" .) The first-order (1PN) expansion of the time component of the local 
equations of motion for a perfect fluid is [12] 

drp + d^) = 0, (5) 



driw + pj+djliw+p + p^ + pu{] = -pd T U, w = p(y + II-[/) (6) 

(U = N.P.[p] is the Newtonian potential associated with the 0-order mass 
density p). The 0-order spatial component of the local equations of motion 
is just the Newtonian equation of motion: 

dripu') + 6>,(p«V) = pUi - p ti . (7) 

Combining the continuity equation (jSJ), and the 0-order expansion of the 
isentropy equation: 

dll = -pd{l/p), (8) 
one gets in a standard way the Newtonian energy equation: 

d T w + dj [(w + p)u j ] = -p d T U. (9) 

Subtracting Q from (0) gives us 

d T p 1 + d J (p 1 u3+pu{) = 0, (10) 

which means that mass is conserved at the first PNA of the scalar theory. 

In contrast, the standard 1PN expansion of the time component of the 
local equations of motion for a perfect fluid in GR consists ([3], Eqs. (1) and 
(64)) of Eq. ©, plus 

d T a + dAou j ) + hpd T U-d T p) = 0, a = p[l + \(u 2 + 2U + U+^)}. (11) 

c A cr p 

(Fock |2] gives only an integral form of the equations of motion.) In the 
equations (0) and (fTTj) of the standard PNA of GR, however, p, u and p 
are not defined except as exact fields ([3], Eqs. (4) and (5)), and, actually, 
they have to be considered as the second approximation of the exact fields — 
instead of being the zero-order expanded fields as in Eqs. (0) and © of the 
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asymptotic PNA of the scalar theory. Thus, each of the matter fields is there 
only in one "exemplar" in the equations of the standard PNA. On the other 
hand, Chandrasekhar |S: does not regard Eqs. © and (|TT|) as exact ones: 
rather, Eq. © is assumed valid if one neglects 0(l/c 2 ) terms, and Eq. (JTTj) 
is assumed valid if one neglects 0(1/ c 4 ) terms. 6 Also note that, unlike 
our Eq. (JHJ) of the asymptotic PNA, the PN equation (jllj) is not a "split" 
equation. Just the same remarks apply to the spatial components of the 
local equations of motion. As to the gravitational field, it is expanded in 
the standard PNA: at the first approximation, there is only the Newtonian 
potential U, whereas "PN potentials" Ui and $ are added in the second 
approximation. The Poisson equation applies to U : 



(Eq. (3) in Ref. j3], Eq. (68.25) in Ref. [2]), and Poisson-like equations are 
derived for [/, and $. 

In our opinion, the difficulty with the standard PNA is this: since the 
matter fields are not expanded, it follows that the equations of the first ap- 
proximation are not exact ones resulting from a "splitting" (as it is the case 
in the asymptotic PNA). Instead, they result from a truncation and can, in- 
deed, be valid only if one neglects 0(l/c 2 ) terms (as Fock and Chandrasekhar 
noted). This point has been proved for the Poisson equation (|12|) in the scalar 
theory, in Ref. [12], Sect. 6.2. Thus, the "asymptotically correct" writing of 
Eq. jnj) of the standard PNA is 



But this means that the field U (the Newtonian potential) can be consid- 
ered to be known only if one neglects 0(1/ c 2 ) terms. Now that field provides 
the main contribution to the acceleration, hence the equation of motion of 
the second approximation (Eq. (68) in Ref. 3) does of course include a term 
involving that field U without an 1/c 2 factor, namely the term pdiU. This 

6 This is clear, e.g., from the sentence after Eq. (109) in Ref. [3] (which is just the 
repetition of Eq. Ijllfl above, Eq. (64) in Ref. 3 ): "That [equation] replaces, in the post- 
Newtonian approximation, the equation of continuity of Newtonian hydrodynamics. We 
may transform the terms in equation (109) that occur explicitly with the factor 1/c 2 with 
the aid of the equations valid in the Newtonian limit." 



AU = —AnGp 



(12) 




(13) 
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means that, at least if one restricts the discussion to the local equations, the 
second approximation does not provide a better approximation than up to 
0(l/c 2 ) terms not included — that is, it does not actually improve over the 
first approximation. Thus, if the first standard PNA is to really reach the 1/c 2 
level included, it must turn out that, at the later stage of the global equations 
of motion for the mass centers, and due to some mysterious integration effect, 
the low accuracy to which the Poisson equation can be considered valid has 
no effect any more. Moreover, as noted by Rendall [5], it may be dangerous 
to have PDE's which are defined only up to unknown higher-order terms, 
which can change the type of the equation, e.g. from hyperbolic to elliptic. 
In conclusion, we prefer to stay with the "asymptotic" PNA summarized 
above, and to see where it leads. 

5 Equations of motion for the mass centers: 
general PN equations 

One may first ask whether there may be any satisfying definition of the mass 
center of a body in a "relativistic" theory of gravitation, given that: (i) in 
a theory accounting for the mass-energy equivalence, any form of material 
energy should be subjected to the action of the gravitational field; (ii) con- 
versely, any form of material energy should contribute to the gravitational 
field; this means that the internal structure (through its energy distribution) 
and the internal motion (through the corresponding kinetic energy) are a 
priori expected to play a role J3j; for these two reasons, it is not obvious 
at all whether one single scalar mass-energy density can be used so as to 
define a single mass center, and which density this should be [15J; (iii) in 
a generally- covariant theory like GR, a covariant definition should be found 
for the mass center — in other words, the motion of the latter should corre- 
spond to a single world-line in space-time, independently of which reference 
frame has been chosen The latter difficulty does not exist in the scalar 
theory, which has a preferred reference frame, but the two first ones subsist. 
Moreover, one has to ensure that the definition adopted will be astronom- 
ically relevant. With modern telescopes and other instruments, the major 
bodies of the solar system are, of course, very-well resolved as extended bod- 
ies, hence the astronomical position refers to an "optical center". It seems 
to be a good approximation to admit that, once corrected from the "phase 
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effects" , the latter corresponds to averaging the rest-mass density, because it 
is the body's "matter", in the usual sense, that does radiate electromagnetic 
energy In any case, this density is certainly more relevant than any den- 
sity involving gravitational energy, for the latter is distributed in the whole 
space. Therefore, we define the (exact) mass center by averaging the (exact) 
rest-mass density in the preferred frame, p ex act 53- A theoretical argument 
also favours that density, and this is the fact that its PN approximation pm , 
at least, obeys the usual continuity equation (see Eqs. (0) and (fTUjl ). thus 
without adding gravitational energy and its flux: this implies that the ve- 
locity of the mass center is itself the barycenter of the velocity — when the 
barycenter is defined with as the weight function [To] . 

Thus, we define the exact mass and mass center through the rest-mass 
density p exact : 

M a exact = I p exact dV, M a exact a cxact ee ! Pcxact xdV (14) 

JD a JD a 

where D a is the (time-dependent) domain occupied by body (a) (a = 1, N), 
in the PRF (V is the Euclidean volume measure in the PRF). At the (first) 
PNA, the mass and the mass center are approximated by 7 

M« = M a + Ml/c 2 , M a ee f pdV, Ml= f Pl dV, (15) 

jD a JD a 

M«a {1) =/ p (1) xrf\/ = M a a + M a 1 a 1 /c 2 , (16) 

JDa 

7 The notation might suggest that, in Eq. 116(1 . the M\a.\jc 2 term is order 2, as the 
product of two first-order terms M\ and ai (still multiplied by the constant coefficient 
1/c 2 ). Recall, however, that the expansions are written in the varying units [M]\ — A[M] 
and [T] a = [T]/-\/A. In these units, M\ and ai are coefficients that do not depend of A, thus 
are order 0, while 1/c 2 is proportional to A. Hence the M^ai/c 2 term is really order 1, not 
2. Similarly, in Eq. (|19|1 . ah terms are order zero. If one wishes, one may come back to the 
starting units (independent of A) , in which the orders of the different fields are not all the 
same, and in which the expansions are hence different, e.g. p cxa ct = A[p + \p\/c 2 ] +0(A 3 ) 
instead of with p and p\ still of order (Ref. Sect. 6). This complicates 

the analysis of the orders, but of course it cannot lead to any inconsistency. Anyway 
in practice, we are then using these expansions for one single value of the parameter A, 
namely the value Ao corresponding to the gravitating system of physical interest (the solar 
system, say). Thus we can use the expressions valid "in the varying units", e.g. (@J. 
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with 

M a a= / pxdV, M 1 a a l = I Pl xdV (17) 

JD a JD a 

Note that M a and a are the Newtonian mass and mass center. Using Eqs. (0) 
and (fTU|). one shows that M a and are constant in time [T5] . To get the 
PN equations of motion of the mass centers, one just integrates the spatial 
components of the PN local equations of motion inside the different bodies 
[T5j . Due to the separation of the different orders in the local equations and 
to their linearity with respect to the PN fields, separate equations are also 
obtained for the mass centers, and the equation for PN corrections (order 1) 
is linear with respect to order-1 quantities. Specifically one finds j!5j : 



M a d l = I pUfdV (18) 

M l a a\ + I ai = J ai + K a \ (19) 
/ [p + p(u 2 /2 + n + U)]uW, (20) 



and 
with 



J m = \ (a^ + pA^dV, (21) 

K ai = / [2k ijPd + P U tl - 2UpUi - a T) k pu'u k - pu'd T K 3 }dV. (22) 

In these equations, a point means (total) derivative with respect to the pre- 
ferred time T; o\ is the PN correction to the active mass density, given by 

a 1 = Pl + p(u 2 /2 + U + U); (23) 

A is the PN gravitational potential, such that the PN expansion of the scalar 
field is 

= 1 - 2U/c 2 - 2A/c A + 0(A 3 ), (24) 

and A is given by: 

A = B + d 2 W/dT 2 , (25) 
where B = N.P.[o"i] is the Newtonian potential associated with o\ and 

W(X,T) ee J GRp(x,T)dV(x)/2 i?=|X-x| (26) 



14 



(G is the constant of gravitation); 2kjj/c 2 and iY l - k /c 2 are the components 
of the non- Euclidean part of the "physical" space metric in the PRF and its 
associated connection, respectively, with 



iFj'fc — Kj,k + Kk,j ^jk,i, Kj — Uhlj, h^- _ ' J ; (27) 



{we use Cartesian coordinates for the Euclidean space metric g°, so that the 
PN physical space metric writes ^2] 



g$ = Sq + } (28) 

and we use Fock's decomposition [2] of any field Z(x), integral of some density 
# x vanishing outside the bodies, into "self " and "external" parts z (x) and 
ZW(x) : 

z(x)= fd x dv = z a (x)+zW(x), ««(x)=/ e x dv,zM(x) = y2 [ e^dv. 

(29) 

The zero-order equation, Eq. (fTHj). is just the Newtonian translational equa- 
tion of motion. 



6 Equations of motion for the mass centers: 
relevant simplifications 

Equation (fT^j) for the PN corrections to the motion is not tractable. To 
make it explicit, we account for two simplifications that do occur for the 
major bodies of the solar system, namely: (i) the good separation between 
bodies (this occurs probably also for many other systems); (ii) the fact that 
the bodies are nearly spherical (this occurs almost certainly also for all other 
systems, if one considers individual bodies with large-enough mass). 



Point (i) is defined by introducing the separation parameter [9\ 

77o = max(r 6 / |a - b|), r b = -Sup x eD |x - y| , (30) 

and by assuming that it is small for the system of physical interest. (The zero- 
order positions of the mass centers are used, since the PN corrections to these 
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positions are very small, hence would lead to nearly the same value for i] .) 
In order to exploit this small parameter, we introduce again a (conceptual) 
family of gravitational systems, by defining initial conditions for them [TU] . 
Thus, remembering the weak-field parameter A, we actually would have a 
two-parameters family of systems. However, the local 1PN equations of the 
asymptotic method — thus the set of the first-order expansions, like (J3J), the 
set of the expanded equations, like (0) and (JTUJ), and the set of the expanded 
boundary conditions, derived in Refs. 12 and 15 — are mathematically self- 
consistent, in that they make a closed system of equations. (They are not 
physically self-consistent yet, since they do not ensure by themselves that 
the weak-field/low- velocity assumptions are satisfied and remain so in time.) 
Hence we may content ourselves with deducing a one-parameter family (S^) 
of well-separated "PN systems" from the data of the given PN system S', 
the latter being itself deduced from the given system S by substituting the 
1PN equations for the exact ones. (Actually the local PN equations of the 
asymptotic method are exact, but they define merely a part of the exact 
fields [12].) The PN gravitational fields: U = N.P.[p], B = N.P.[<7i], and W, 
depend merely on the PN matter fields [Eqs. I2*1H2*o] , Hence we just have to 
define the initial conditions for the PN matter fields. Moreover, the latter 
conditions turn out to be determined by the initial zero-order matter fields 
p(T = 0), or equivalently p(T = 0), and u(T = 0) [HJ. To define p(T = 0) 
in S^, we set 

a"(T = 0) = a(T = 0) Vo / V , (31) 

p 7? (x,T = 0) = p(a + y,T = 0) if x = a" + y with a + yeD a . (32) 

Equation defines p v so that the density inside the bodies is independent 
of r) [setting p J? (x, T = 0) = if x does not have the form above for some 
a = 1, N]. In other words, the bodies themselves do not depend on the 
separation parameter 77. Equation (|31|) ensures that, at least near T = 0, the 
separation distances between bodies are of order rf 1 : 

| a "-b ,? | =oid{r ] - 1 ). (33) 

To define the velocity u v (T = 0), we use the auxiliary assumption that each 
body undergoes a rigid motion at the Newtonian approximation [2]: 

u l = a' + {x ] — a/), or u = a + uj a A (x — a) for x G D a , + fl\f = 0. 

(34) 
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Of course, this is only approximately true (e.g. due to the tidal influence of 
the other bodies), but we use this assumption merely to calculate the PN 
corrections. Since the latter ones are very small, it certainly implies only 
an extremely small error in the solar system. Anyhow, we can assume that 
is exact at the initial time. From the Newtonian estimate a 2 « [/^(a) 
, valid in the reference frame of the global mass center, and from (J33)) . we 
expect that, at any time, 

(a*)" = ord(?? 1/2 ). (35) 

We assume that this is true if we define the initial translation velocities of 
system S v as 

( a y(T = 0) = (r ] /r ]0 ) 1 / 2 d l (T = 0). (36) 

As to the self-rotation velocities, in the solar system they have at most the 
same magnitude, in linear values, as the translation velocities, and our nu- 
merical calculations show that the PN corrections containing quadratic terms 
in ilji , included in the first version of the explicit translational equations [§] , 
are negligibly small in the solar system. To avoid such terms in the expan- 
sions, it turns out to be sufficient that 

(^ ) )"«V /2 , (37) 
hence we set, for some small number e > 0, 

(nW)"(T = o) = ( v /vo) £+1/2 n«(r = o), (38) 

and we assume that this ensures that (|3*7)l is true at any time. 

Point (ii) is defined simply by assuming, merely at the stage of calculating 
the PN corrections, that the zero-order rest-mass density p is spherically 
symmetric for each body: 

Vx e D a , p(x) = p a (r), r=|x-a| (a — 1, N). (39) 

The sphericity of the field p implies also that of the pressure field p and the 
Newtonian self-potential u a . By using this and point (i), i.e. Eqs. (JHHJ), (|3*3jl 
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and (|37|) , we get after rather involved calculations [TDJ : 



du la 
dT 



with 



o(v 3 ) + 



\(- 









M, 



^ + 2^ u Q .VC/W(x a 



17e a + 6T a 
3M a 



2— + /7( a Vx 



V[/ (a) (Xa) 



+G y^ / Mb_ [( na6 _u 6 ) na6 _ 



b^a 

+ — i -®ab n ab + ~ 

T ab 

+ [xi 6 -xi a + 3 ((xi a -x 16 ).n a6 ) n ab ] ) , 

r ab 

|x a - x 6 |, n a6 = (x a - x b )/r ab , and 



(3(n ab .u 6 ) 2 - uf) n afe - (n afe .u fe ) ( 2u fe + -u c 



(40) 



"16 = M b 



uM(x a ) + U®{x b ) + 



ll£a + 8T a 
3M„. 



(41) 

and where, as a preparation for the application, the notation for the PN 
positions and velocities of the mass centers has been changed, as compared 
with (|TH|I and (fT7 |l : 



x a = a, Xi„ = c (a ( i) 



(l)o = a (l) = X « + Xla/C 2 (42) 



u a = a, uia = c 2 (a (1) - a) = x lo , u ( i )o = a (1) = u a + u la /c 2 (43) 

Note that terms of order up to and including rj 3 have been consistently re- 
tained in Eq. (|4()jl [10J. In the first version of the explicit translational equa- 
tions of motion [§], this could not be done due to the lack of an asymptotic 
framework for the separation parameter rj, and this resulted in a large dif- 
ference with observational data. In addition to the self-rotational energy: 

T a = 4° = I P(x* - a')(xi - a^dV, (44) 

JD a 

two structure- dependent parameters appear in ()4(Jj) : 

e a = f pu a dV/2, (45) 

JD a 
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3G 



| dr, n a (r) = 4tt J p a (s)s 2 ds. 



(46) 

No such structure parameter does enter in the PN equations that have been 
used in celestial-mechanical tests of GR, namely the Lorentz-Droste-Einstein- 
Infeld-Hoffmann (LDEIH) equations, which are based on the standard PNA. 
However, as already noted by Synge and coworkers, one should expect that 
the internal structure of the bodies does influence the motion. 



7 Numerical implementation. Comparison with 
an ephemeris based on the standard PNA 
of GR 

In order to use the translational equations of motion (fTHj) and (|4*0|) . so as to 
check the scalar theory, we have to know the values of the independent pa- 
rameters that enter these equations. These are: the zero-order masses M a of 
the bodies (here the major bodies of the solar system) and their parameters 
T a , e a and r a ; the initial conditions of their motion; and the constant velocity 
V of the global zero-order mass center of the solar system, with respect to 
the preferred frame E (Of course there is no parameter like V in con- 
ventional theories.) These unknown parameters depend on the theory. They 
have to be determined by optimizing the agreement between predictions and 
observations [9]. In particular, one may expect that the optimal values of 
the zero-order parameters should differ from their values in pure Newtonian 
theory by first-order (second-approximation) quantities, this fact alone giv- 
ing corrections of the same magnitude as the 1PN corrections jH]- However, 
since the parameters T a , e a and r a play a role only in the PN corrections, 
such first-order corrections on their values would make only second-order dif- 
ferences, hence negligible ones, in the final results. Therefore, we calculate 
these parameters from the standard rotation velocities and density profiles 
of the corresponding bodies j2H • 

Then our adjustment algorithm loops on the numerical solution of the 
translational equations of motion in order to optimize the remaining parame- 
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ters, i.e., the zero-order masses M a , the initial conditions, and the velocity V. 
The algorithm has been tested by investigating in which measure one may 
reproduce (over one century) the predictions of the DE403 ephemeris j^Z] , by 
using purely Newtonian equations of motion. (The DE403 ephemeris is based 
on LDEIH-type equations of motion [28J, thus it is based on the standard 
PNA of GR.) Our algorithm has also been applied to adjust a less simplified 
model, in which the PN corrections in the Schwarzschild field of the Sun are 
also considered, and to compare this model with DE406 over 60 centuries 
30J. It has thus been found that the difference between a calculation based 
on these "Schwarzschild-corrected" Newtonian equations, limited to the Sun 
and the nine major planets, and the DE406 ephemeris of the JPL, based on 
LDEIH-type equations and including the Moon and asteroids, is very small. 
Over the last century, for instance, the longitude difference for Mercury is 
0.03" [20]. The VSOP82 ephemeris IH] is also based on "Schwarzschild- 
corrected" Newtonian equations, but it uses more accurate (semi-analytical) 
algorithms than ours; also, for the comparison with the JPL ephemerides, it 
seems that it takes Schwarzschild's metric in "isotropic" form, rather than 
in the standard form in which we took it. The VSOP82 ephemeris leads to 
even smaller differences with an ephemeris based on LDEIH-type equations 
(e.g. 0.001" over 1891-2000 for the longitude of Mercury). Now the scalar 
theory does predict Schwarzschild's motion for test particles in the field of 
one spherical body, if the latter is at rest in the preferred frame (V = 0). 
Hence, if we take for PN corrections in our theory those that are obtained by 
considering the planets as test particles in the field of the spherical Sun, and 
if we include V in the free parameters, then we can obtain only a smaller 
difference with the LDEIH-type equations than the difference between the 
latter and the "Schwarzschild-corrected" Newtonian equations. But, in our 
opinion, the correct equations of motion are those got with the asymptotic 
PNA. Therefore, we have implemented the equations ([18)1 and (|4U|) in our 
adjustment algorithm. 

To do that, we had first to solve a numerical shortcoming of the asymp- 
totic method, namely the fact that the mass center of a body, or the test 
particle, is followed on its trajectory by two different positions: the zero- 
order position xo and the 1PN position, X(i) = xo + xi/c 2 , which drift from 
one another as the time goes. By investigating in detail the case of a test 
particle in a Schwarzschild field, it has been found that this leads to an error 
increase in (T — T ) 2 , and that it may be cured by "reinitializing", i.e., by 
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substituting T + 5T, then T + 2ST, etc., with \8T\ sufficiently small, for the 
initial time T Q in the differential system that governs the motion of the mass 
centers [32j • Moreover, since the equation for PN corrections (J4*U|) is valid 
only in the preferred reference frame E, we use a Lorentz transform to go 
from E to the frame Ev bound with the zero-order global barycenter, and 
vice- versa. This transform, as well as the inverse transform, is determined by 
the free vector V. Thus, the adjustment process of the translational equa- 
tions on observational data provides us eventually with the value of V that 
minimizes the residual with the set of observations. However, the "observa- 
tional data" are currently taken from an ephemeris based on GR, specifically 
here we took from DE403 a set of heliocentric positions of the eight major 
planets (Pluto omitted), between 1956 and 2000. 

With these input data, themselves a fitting of observations by the LDEIH 
equations of GR, the magnitude of the optimal vector is |V| ~ 3km/s, which 
is significant. The difference between DE403 and our thus-adjusted equations 
of motion is shown on figures 1 and 2. The self-rotation of all nine bodies is 
neglected, i.e., all uj a 's are assumed in (J34j) . We show only the difference 
in longitudes, because it dominates over the other errors. It can be seen 
that, for most planets, the difference is very small (a few times 0.1") over 
the fitting interval, but it increases quickly with time for the inner planets 
(Mercury, Venus, the Earth- Moon Barycenter (EMB), and Mars). We do 
not know yet whether this comes only from the different models or partly 
also from numerical reasons: due to the necessity of reinitializing very often 
(every two days here — we mean the ephemeris time, not the computer time), 
the calculations are long. It has already been checked that an increase of the 
accuracy in the ODE-solver brings negligible changes to the present Figure 
showing the differences between DE403 and the asymptotic PNA of the scalar 
theory. Anyhow, the differences are still quite small, e.g. 3.8" for Mercury 
after the last century, to be compared with the relativistic perihelion advance 
of 43", and with the accuracy of the current ephemerides, considered to be 
0.1" for Mercury [HII|,p. 228. For the influence of the Moon on the motion of 
the EMB, we use a semi-analytical correction formula j33], which we adapted 
in an approximate way from the "1950 ecliptic" to the "J2000" reference. 
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Figure 1: Longitude differences (arc seconds) for planets Mercury to March, ei- 
ther as taken from the DE403 ephemeris of the Jet Propulsion Laboratory, or as 
obtained (after adjustment) by numerical integration of the equations of motion in 
the asymptotic post-Newtonian approximation of the investigated scalar theory. 
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Figure 2: Longitude differences (arc seconds) for planets Jupiter to Neptune, 
either as taken from the DE403 ephemeris of the Jet Propulsion Laboratory, or as 
obtained (after adjustment) by numerical integration of the equations of motion in 
the asymptotic post-Newtonian approximation of the investigated scalar theory. 
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8 Conclusion 



The usual method of asymptotic expansions, as defined in Sect. 2, is in- 
deed of usual utilization in most domains where partial differential equations 
occur, but it contrasts with the standard (Fock-Chandrasekhar) method of 
post-Newtonian approximation (PNA) for weak gravitational fields. In the 
latter method, no one-parameter family of similar problems is introduced, 
so that the meaning of the approximation is not very clear. Indeed 1/c 2 
is formally considered as a small parameter, and the matter fields are not 
expanded. We applied the usual method of asymptotic expansions to weak 
gravitational fields in our scalar theory, and we call the result the asymptotic 
PNA. A similar method has been proposed in GR by Synge and coworkers 
[SJ 0] for stationary gravitational fields, and has been initiated in a more gen- 
eral case, but not fully developed, by Futamase & Schutz jTj. In our opinion, 
the asymptotic PNA is very solid mathematically. Within the asymptotic 
PNA, 1/c 2 turns out to be (proportional to) the small parameter A, but 
this is true in specific units, depending on A. The asymptotic PNA leads to 
definitely different equations, as compared with the standard PNA. In par- 
ticular, the former method predicts that the internal structure of the bodies, 
and their internal motion, has a definite influence on the motion of the mass 
centers of a self-gravitating system of bodies. 

This has been checked numerically in the solar system for the scalar the- 
ory. The standard PNA should probably lead to the same result in the scalar 
theory as in GR, namely it should lead to the conclusion that, in the solar 
system, the post-Newtonian effects may be calculated simply by adding to 
the Newtonian motion the PN corrections obtained in considering all planets 
as test particles in the field of the isolated Sun. If the latter approximation 
is used, the results of our scalar theory are nearly indistinguishable from 
those of GR. 8 In contrast, if the equations of the asymptotic PNA are used, 
then the predicted motion apparently cannot fit a standard ephemeris (i.e. 

8 We did check this numerically. To this end, for each planet, we took the PN equation 
of motion of a test particle |35) (p. 22), the first-order contribution assuming a spherical 
Sun. The equation of motion then reduces to the "Schwarzschild-corrected Newtonian 
equation" , plus extra terms which vanish if the velocity of the Sun through the preferred 
frame is zero. As a result of fitting these equations to the DE403 ephemeris, that velocity 
was found negligible, and of course its effect on the thus-calculated ephemeris was then 
found negligible also. 
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an ephemeris based on the standard PNA of GR) within what is currently 
believed to be the observational accuracy. 

It seems that the influence of the internal structure of the bodies and 
the difference with the DE403 ephemeris are much more the result of chang- 
ing the approximation method, than that of changing the theory. Indeed, 
the exact local equations of motion for a perfect fluid in the scalar theory 
are very similar to those in GR, and the general metric of the theory is a 
"Schwarzschild-like" metric. In other words, the author considers it likely 
that a similar departure from a standard ephemeris would be left, if one 
compared it with a calculation based on an asymptotic PNA of GR. How- 
ever, this could be proved only by building a general asymptotic scheme in 
GR. This would be difficult, due to the fact that in GR the initial conditions 
cannot be imposed freely. 

Coming back to the scalar theory, there is still a possibility to improve the 
fitting, mainly by re-adjusting the masses (this could not be done here, except 
for those of the Sun and Jupiter, because the masses are sensible parameters 
whose values have to be separately optimized before possibly optimizing them 
globally), perhaps also by improving the numerical accuracy. It is also very 
important to adjust the equations, not on an ephemeris (because this is 
already a fitting of observations by some other equations), but directly on 
observations. Indeed some correction factors of observational data are taken 
as free parameters in the adjustment of an ephemeris, hence the observations 
are not completely independent of the gravitational model. Finally, it should 
be noted that the (best-fitting) value of the absolute velocity V of the mass 
center of the solar system has been found to be ca. 3 km/s with the present 
model (self-rotation neglected, adjustment on a standard ephemeris), and 
this already is not negligible. Due to these simplifications, and due to the 
fact that the solar system was assumed isolated here, this present best-fitting 
value of V is not even an approximation to the correct one: it justs tells an 
idea about the order of magnitude of V. That this is not negligible, might 
incline one to think that the preferred-frame character of the theory is not 
redhibitory. 
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